Slow dynamics of a confined supercooled binary mixture: direct space analysis 
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Dynamical properties of a Lennard- Jones binary mixture embedded in an off lattice matrix of 
soft spheres are studied in the direct space upon supercooling by molecular dynamics simulations. 
On lowering temperature the smaller particles tend to avoid the soft sphere interfaces and corre- 
spondingly their mobility decreases below the one of the larger particles. The system displays a 
dynamic behaviour consistent with the Mode Coupling predictions. A decrease of the mode cou- 
pling crossover temperature with respect to the bulk is found. We find however that the range of 
validity of the theory shrinks with respect to the bulk. This is due to the change in the smaller 
particle mobility and to a substantial enhancement of hopping processes well above the cross over 
temperature upon confinement. 

PACS numbers: 61.20.Ja, 61.20.Lc, 64.70.Pf 



I. INTRODUCTION 

Liquids under confinement represent a field of grow- 
ing interest in science due to the connection with rel- 
evant technological and biophysical problems^. In this 
field modifications of both the phase diagram and dy- 
namics of the confined liquid with respect to the bulk 
represent a key point, in particular as far as the possibil- 
ity of supercooling is concerned?. Experiments exhibit a 
very diversified phenomenology for supercooled confined 
liquids. In particular the calorimetric glass transition 
temperature, T g , can increase^, decrease^ or be unaf- 
fected^ with respect to the the bulk value depending on 
the liquid, the confining geometry and the nature of the 
substrate. 

Phcnomenological arguments predict the existence of 
domains of cooperative dynamics in the supercooled liq- 
uidi. The size of these domains is expected to grow upon 
supercooling. The existence of an upper bound for the 
domain size in confined media implies a decrease of T g 
upon confinement. In fact if this size remains small on 
supercooling then the probability of a cooperative mo- 
tion, and therefore of a structural relaxation, is large, 
and dynamics is faster. 

The mode coupling theory of the glass transition 
(MCT) 8 is able to describe the dynamics of bulk liquids in 
the supercooled region on approaching a crossover tem- 
perature Tq- At Tc the system passes from a regime 
where ergodicity is attained through structural relax- 
ations to a regime where this mechanism is frozen and 
only activated processes permit the exploration of the 
configurational space. Tc can be estimated through ex- 
periments and computer simulations or predicted by the 
theory. 
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Above Tc the relaxation mechanism of the supercooled 
liquid can be described as mastered by the cage effect. 
Nearest neighbors surround and trap the tagged parti- 
cle forming a cage around it. When the cage relaxes, 
due to cooperative motions, the particle moves. MCT 
translates this phcnomenological description into a pre- 
cise mathematical framework where, starting from the 
classical equation of motion for the density correlator 
4>q{t), a retarded memory function is introduced. In the 
idealized version of MCT hopping processes are neglected 
and the retarded non linear set of integro-diffcrential 
equations can be solved analytically to the leading or- 
der in e = (T — Tc)/Tq, deriving universal results for the 
density correlator. With these approximations Tc is the 
temperature of structural arrest of the ideal system. The 
success of this theory is due to the fact that on approach- 
ing Tc from above hopping processes can be neglected in 
many liquids and the predictions of the idealized version 
of MCT hold. 

Molecular dynamics studies of model liquids in re- 
stricted geometries intended to assess the applicability of 
MCT can give an important contribution to the charac- 
terization of the glass transition scenario in confinement. 
In this field recent progresses have been done for water 
confined in a silica nano por o 9 i 10 and for confined thin 
polymer films 11 . 

We performed a Molecular Dynamic (MD) study of 
the single particle dynamical properties of a Lennard- 
Jones binary mixture (LJBM) embedded in an off lattice 
matrix of soft spheres. The model considered represents 
a situation of strong confinement analogous to real silica 
xerogelsi2iiii4. The bulk phase of the chosen mixture 
is known to test most of the MCT predicted features^. 
A brief report on some of our findings has been recently 
published^. 

We present in this paper direct space quantities eval- 
uated from MD trajectories in order to define the range 
of validity of MCT with respect to the bulk phase. In 
the next section we describe the model adopted and the 
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details of the MD simulations. In the third section we 
illustrate the radial pair distribution functions and the 
coordination numbers of the system. The fourth section 
is devoted to the study of the dynamic quantities in the 
direct space: the mean square displacement (MSD), the 
van Hove self correlation function and the non-Gaussian 
parameter. Last section contains the concluding remarks. 




FIG. 1: Graphic representation of the model studied. Small 
light gray particles are of A type, small dark gray particles are 
of B type. The soft spheres that form the confining system 
appear in the figure as large dark gray spheres. The cube 
in which the system is embedded is the simulation box. The 
radii of the three types of particles are purely representative. 
The a values are listed in TabQ] 



II. MODEL AND MOLECULAR DYNAMICS 

The confining medium consists in a rigid disordered 
array of soft spheres. The simulation box contains 16 
soft spheres labeled in the following with the letter M. 
The interspaces of this structure host a liquid LJBM of 
1000 particles, composed of 800 particles of type A and 
200 particles of type B. The interaction pair potential 
can be written in the following general form: 



V lu ,{r)=4e„ 



(i) 



where indexes /i, v run on the particle types A,B,M. The 
parameters for the three components are listed in Ta- 
ble 0] The parameters of the LJBM have been chosen 
as in refpi^. The shifted potential technique has been 
used with a cutoff of 2.5 a^ v . Periodic boundary condi- 
tions are applied. In the following Lennard- Jones units 
will be used, therefore energy will be given in units of 




FIG. 2: Graphic representation of the confining disordered 
structure. The sixteen spheres of the matrix contained in the 
simulation box of FiglTlhave been repeated to obtain a box of 
height 2L, width 2L and depth L. The surface that appears 
in the figure is equipotential. 



temperature in units of 6aa/Kb, length in units 
oaa and time in units of (m(j\ A / (ASseaa)) 1 ' 2 ■ The box 
length is L = 12.6. The simulation box of the system 
is represented in Fig. ^ As it can be seen from the pic- 
ture the LJBM is in a strongly confined environment as 
only few layers of particles can accommodate among the 
soft spheres. Fig. [21 shows a cell made of four simulation 
boxes. Only equipotential surfaces of the soft spheres are 
displayed in order to offer a clearer image of the disor- 
dered porous structure hosting the liquid and to highlight 
the analogies with silica xerogels. 

We have conducted MD simulations of this system in 
the NVE ensemble. The equations of motion were solved 
by the velocity Verlet algorithm. The system was equi- 
librated at different reduced temperatures via a velocity 
rescaling procedure. The temperatures investigated are 
T = 5.0,2.0,0.8,0.58,0.538,0.48,0.465,0.43,0.41,0.39. 
The timestep used for T > 1 was 0.01 and for T < 1 was 
0.02. For the lowest temperature investigated, T = 0.37, 
a production run of 14 x 10 6 timesteps was performed. 
We verified that the results obtained do not depend on 
the specific choice of the disordered matrix by running 
MD simulations for two other different configurations of 
the disordered matrix of soft spheres. Both thermody- 
namics and dynamics appeared the same for the three 
systems, in particular deviations among the diffusion co- 
efficients extracted at given temperatures from the MSD 
of different configurations are within 2%. 

Values of thermodynamics quantities of the equili- 
brated systems are shown in Fig|3J Averaged pressure, 
potential energy and total energy per particle of the mix- 
ture behave similar to the bulk. A stability of the system 
throughout all the isochore explored is evident from the 
smoothness of the curves. 
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Due to the soft spheres potential the free volume acces- 
sible for the mixture strongly depends on temperature. 
In Fig. 01 we show an estimate of the free volume for A 
and B particles in our confined system as a function of 
temperature. The estimate has been carried out through 
the Voronoi tessellation. This procedure has been ap- 
plied only to non-interfacial particles. The free volume 
has been calculated by averaging over 160 configurations 
of the system for each temperature. 
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FIG. 3: Thermodynamics of the confined Lennard- Jones bi- 
nary mixture: pressure (circles), total energy per particle 
(squares) and potential energy per particle (diamonds) ver- 
sus the inverse of the temperature. 
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III. RADIAL PAIR CORRELATION 
FUNCTIONS 

The radial pair correlation functions g^v{r) are shown 
in Fig|3J For all the temperatures investigated the system 
is in the liquid phase as no abrupt sharpening of the peaks 
is evident. No indications of phase separation appear 
upon supercooling either. The latter phenomenon would 
in fact result in a substantial decrease of the area of the 
first peak, ab\ , of the g^B (f) , which we do not observe in 
our system. On lowering temperature the peaks sharpen 
and enhance more than in the bulk due to the strong 
packing induced by the soft spheres. A clear element of 
distinction with respect to the bulk is the enhancement 
of intensity of the first gBB(r) peak, bbi, as supercooling 
progresses. The same peak appeared instead reduced to a 
small shoulder for supercooled temperatures in the bulk. 

The functions guA (f) and gm b (f) are characteristic of 
this confined system and are representative of the interac- 
tion between the confining soft spheres and the particles 
of the mixture. We note the gradual disappearance of 
the first peak, mbi, of the gMB(r) that begins at T — 2.0 
while the second peak, 77162, enlarges its area. Corre- 
spondingly the first peak of the gMA{r) shifts to larger 
r. For lower temperatures, namely T < 0.8, the struc- 
ture of the mixture stabilizes so that MA and MB peaks 
alternate in space. 

In order to quantify the above considerations for the 
radial distribution functions we calculated the coordina- 
tion number defined as: 



(2) 



where p is the total number density N/V, where N 
1 01 (i. V — L 3 and x v is the relative fraction of the v 
species. The quantity of Eq. [21 evaluated for each peak 
of the g^v (r), is shown in Fig. [5] The number of B parti- 
cles around the soft sphere in the first shell lowers upon 
decreasing temperature in favor of an enhancement in the 
second shell. Correspondingly the number of A particles 
in the first shell enhances. 

The above considerations on the radial distribution 
functions and the coordination numbers lead to the con- 
clusion that B particles avoid the interface with the soft 
spheres for low temperatures. In Fig. [7| a snapshot of 
the system for T = 0.37 evidences this behavior. The 
fluctuations of the error bars in Fig0| also confirm this 
phenomenon. In fact since interfacial particles have been 
excluded from the calculation of the free volume a larger 
fluctuation corresponds to a frequent exchange of the par- 
ticles with the interface. 



FIG. 4: Specific volumes of the Voronoi cells for A and B 
particles. Particles in contact with the soft confining spheres 
have been excluded. The bulk considered is that of ref^. 



IV. DYNAMIC QUANTITIES 

In the following we shall describe the temperature vari- 
ations of the time dependent direct space quantities ana- 
lyzed for our confined system upon supercooling, namely 
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FIG. 5: Radial distribution functions calculated at differ- 
ent temperatures for the system. From top to bottom AA, 
BB, AB, MA and MB. Lower temperatures are progressively 
shifted upward. For the AA function the shift is 0.2, for AB 
0.3 while for the remaining BB, MA and MB it is 0.1. The 
labeling of the peaks is used in the text and in Fie'151 

the MSD, the van Hove self correlation function and the 
non-Gaussian parameter. 

The asymptotic behaviour of the MSD can be used as a 
test of an MCT prediction. The MSD is in fact expected, 
as supercooling progresses, to develop a plateau at inter- 
mediate time scale, corresponding to the rattling of the 
particle in the cage. After the initial short time ballistic 
motion the particle enters the plateau whose extension 
in time enhances upon lowering temperature. After leav- 
ing the plateau region the Brownian behaviour predicted 
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FIG. 6: Coordination numbers of the shells defined by the 
peaks of the radial distribution functions as labeled in Figl5l 
See Eq. H 




FIG. 7: Tomographic snapshots of the system at the lowest 
temperature investigated, T=0.37. In order to obtain this pic- 
ture the system has been divided in slices of height 1.4a. From 
top to bottom and from left to right every picture displays 
the content of each slice. Empty circles are A-type particles, 
filled circles are B-type. Crosses mark the volume of the soft 
spheres. No phase separation is evident. At this temperature 
B particles avoid the interfacial region. 

for a simple liquid in a stable state is restored. From the 
slope of the late MSD the diffusion constant D can be ex- 
tracted through Einstein relation < r 2 (t) >= 6Dt. The 
temperature behaviour of diffusion coefficients can give 
an estimate of the mode coupling crossover temperature 
through the MCT power law behaviour 

D ~ (T - T c y (3) 

In Fig|H]we show the MSD calculated for both A and B 
particles as a function of time for all the temperatures 
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FIG. 8: Log-log plot of the mean square dis- 
placements of A particles and B particles for 
all the temperatures investigated, namely T = 
5.0, 2.0, 0.8, 0.58, 0.538, 0.48, 0.465, 0.43, 0.41, 0.39. 



investigated. It can be seen that also in this confined 
system the binary mixture displays the typical features 
of a glass former. From the height of the plateau the cage 
radius can be extracted and the value obtained for both 
A and B particles is r 2 = 0.04, analogous to the bulk 
result. 

The lowest temperature reached in our study is T = 
0.370 against T = 0.466 for the bulk. The elongation 
of the plateau of the MSD for our lowest temperature 
is comparable to that of the lowest temperature of the 
bulk. The diffusion coefficients extracted from the slopes 
are reported in Fig. At higher temperatures T > 1.0 
the larger A particles diffuse slower than B particles while 
they move faster at lower temperatures. As shown in the 
previous section, in the confined case for lower temper- 
atures B particles tend to avoid the interface of the soft 
spheres. As a consequence they remain confined in the 
inner part of the channels, surrounded by A particles, 
and this diminishes the accessible free volume lowering 
mobility below the one of A particles, see also Fig0J 

In Fig. the fit to the MCT power law of Eq. is 
also reported. The estimate of the crossover temper- 
ature given by the fit is Tc = 0.356 for both A and 
B particles. The exponents are 7 = 1.86 and 1.89 re- 
spectively for A and B particles. We do observe devi- 
ation from a power law both at high and low tempera- 
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FIG. 9: Diffusion coefficients for A and B particles (symbols). 
Continuous lines correspond to the fit to Eq. [3] The values 
of the parameters extracted from the fit are Tc = 0.356 for 
both species and 7 = 1.86 for A particles and 7 = 1.89 for B 
particles. 



tures. The temperatures used for the fit are in the range 
0.410 < T < 0.58 corresponding to 0.153 < e < 0.631. 
The range of validity of MCT reported for the bulk is 
much wider, 0.07 < e < 1.30. The upper limit is gen- 
erally due in bulk to the non-validity of the asymptotic 
expansion for too large values of e. Here we infer that a 
further decrease of this limit is due to the progressive B 
particles avoidance of the soft sphere interfaces upon su- 
percooling. For the temperatures below the upper limit 
the peaks of g(r) for MA and MB have in fact already 
stabilized, see FigO The deviation observed at the two 
lowest temperatures is to be connected to the presence of 
hopping mechanisms, analogously to bulk liquids, as we 
will see in more detail when we illustrate the behaviour 
of the VHSCF. 

The Van Hove self correlation function, VHSCF, is de- 
fined, for a system of N particles, as 



1 N 

G s (r,t) = -(£S[T + r i {0)-T i (t)]) 



(4) 



4irr 2 Gs(r,t)dr is the probability of finding a particle at 
distance r after a time t if the same particle was at the 
origin r — at the initial time t = 0. 

In Fig|n3 we show 4Trr 2 Gs(r,t)dr at fixed times for 
three different temperatures. On this scale the correla- 
tors of A and B particles do not differ substantially. At 
high temperature the VHSCF decays regularly as time 
increases as expected for a simple liquid. As supercool- 
ing progresses the system begins to display a new feature 
that is best evident for the lowest temperature. This 
feature consists in the clustering of the curves at inter- 
mediate times which is the signature of the caging of 
the particles. The zone where the peak value clusters is 
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r = 0.15 — 0.25 and corresponds to the zone in which the 
MSD flattens. For the lowest temperature investigated 
the maximum and minimum times for the clustering are 
2.6 < t < 82.0. These values define the so called (3 re- 
laxation temporal region. In the comparison between 
the VHSCF of A and B particles it can be observed the 
crossing between the velocity of the two species. At high 
temperatures, in fact, the VHSCF of B particles extends 
further than the one of A particles for long times, while 
for low temperatures behaviors are swapped. 




FIG. 10: 4-Kr 2 Gs(r,t)dr vs r. a) b) and c) pictures refer to 
particle A for T=5.0, 0.5 and 0.48 respectively. Data are sam- 
pled at fixed times. The times chosen follow the progression 
t — 2". In the pictures d) e) and f) the correlators are shown 
for B particles (dashed lines) together with the correlators for 
A type for a comparison (continuous line) . The curves chosen 
for a comparison are those indicated with a thicker line in the 
graphs on the left. 

A behavior not present in the bulk can be observed 
for the lower temperatures by looking at the peak of the 
VHSCF for times longer than the plateau region. In this 
temporal range the peak of the VHSCF is expected, in 
the framework of MCT, to shift to larger r for large times. 
The absence of this shift in our confined mixture, as we 
shall see in the next picture, must be connected with an 
important hopping phenomenon able to restore diffusive 
motions when the cage is not jet dissolved. 

In Fig^Jwe show a blow up of the tails of the curves 
of Fig^H Aside the previously cited cage peak, we note 
the presence of a hopping peak around r — 1.0. This 
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FIG. 11: Blow up of the 4-Kr 2 G s {r, t)dr vs r for T=0.37 for 
different times for A (top graph) and B (bottom graph) parti- 
cles. Evidence is put on the region in which the peaks related 
to the presence of hopping appear. 



peak enhances for late time and it is more pronounced 
for A particles than for B. In the bulk on approaching Tc 
from above a hopping peak appears only for B particles. 
From the MSD it is evident that our hopping phenom- 
ena can not be neglected at the two lowest temperatures, 
namely at T = 0.37,0.39. As a consequence these tem- 
peratures can not be accounted for a test of the MCT in 
the region around the cage relaxation time, a relaxation 
region. Nevertheless these temperatures can be used to 
test the theory in the /3-time region. In fact, see also 
ref^, hopping peaks do not affect curves in the time 
range of the clustering region. 

For the A particle a second broader peak is visible in 
Fig^2 This peak is also connected to the presence of 
hopping. The two hopping peaks correspond to the first 
and second peak of the gAA(f), see FiglSJ meaning that 
there are some privileged position for the atoms to be 
reached and that those atoms that reached these posi- 
tions have traveled more than the average. Similar, the 
hopping peak of the B VHSCF corresponds to the first 
peak of the g_Bs(r). Hopping phenomena are connected 
with the existence of dynamical heterogeneities in super- 
cooled liquidsii. In confinement they appear to have 
above Tc a more important role than in the bulk. 

Dynamical heterogeneities are also known to cause a 
non gaussianity of atomic displacements. We evaluated 
the first non-Gaussian parameter (NGP), that works as 
a quantifier of the degree of non-Gaussianity. This pa- 
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FIG. 12: Non Gaussian parameter calculated for A (top) and 
B (bottom) particles for all the temperatures investigated. In 
both graphs curves on the top correspond to lower tempera- 
tures. 



rameter is defined^ in three dimension as 
3 < r 4 (t) > 



a 2 (t) 



5 < r 2 (t) > 2 



1 



(5) 



This quantity calculated for our system is reported in 
fig El Also for this parameter we have marked differ- 
ences with respect to the bulk that reflect the different 
weight of the hopping processes in the two cases. In the 
bulk the NGP of A and B particles were different and 
never exceeded the values of 2 for A particles (no hop- 
ping above Tc) and 3.2 for B particles. Here instead we 
have similar NGP for A and B particles and the value 
grows enormously as we supercool. For the last tem- 
perature investigated the peak is around 5. We found 
no relation between the peak positions and MCT power 
law of EqO as it was found for example in bulk water— 
and the curves of the confined NGP do not appear to 
follow a master curve before the maximum as found in 
bulk liquidsi^ia. The trend of the NGP of bulk and 
confined mixture depicts a scenario in which MCT be- 
haviour causes a mild deviation from gaussianity. The 
peak is around 2 when close to Tc^^, its position shifts 
substantially upon supercooling and the NGP follows a 
master curve in the temporal region before the peaks. 
Hopping appears to cause a marked deviation from gaus- 
sianity instead, leading to a disappearance both of the 
MCT master curve and of the shift of the peaks upon 



supercooling. 

V. SUMMARY AND CONCLUSIONS 

We presented an MD study of a LJBM confined in a 
disordered array of soft spheres. The direct space quan- 
tities have been analyzed upon supercooling for a test of 
the MCT behaviour. 

Two main differences with respect to the bulk mixture 
appear due to confinement. The g(r) have evidenced 
that smaller B particles tend to avoid interfaces at low T 
and correspondingly their diffusion coefficient becomes 
lower than that of the A particles. Hopping processes 
are present also for A particle above Tc as shown in the 
VHSCF. 

Not obviously, these differences do not prevent the 
MCT asymptotic predictions from holding for this mix- 
ture also in this kind of confinement. The range of valid- 
ity of MCT in the region of the a relaxation is how- 
ever much more limited. The range of the bulki£ is 
0.07 < e < 1.30 against the 0.153 < e < 0.631 found 
in our confined model. Disturbances in the dynamical 
behaviour of the bulk introduced by a different type of 
confinement are also reported to lead to a complete dis- 
appearance of MCT—. 

Power law fit to the diffusion coefficient extracted from 
the slope of the MSD give an estimate of Tc = 0.356 for 
both species and 7 = 1.86 for A particles and 7 = 1.89 
for B particles in agreement with MCT that states that 
these quantities should be the same for both A and B 
particles. The values of the bulk extracted from D in 
literature are±£ T c = 0.435 for both species and 7 = 2.0 
for A particles and 7 = 1.7 for B particles. The exponents 
7 appear similar while we have a reduction of Tc upon 
confinement. A complete test of the correlators in the Q 
space will be presented in a subsequent paper™. 

Acceleration of dynamics in presence of repulsive con- 
finement is expectedii based on the fact that in confine- 
ment the cooperative rearranging region cannot grow be- 
yond the size of the system. In the case of an attractive 
interaction this effect might superpose to the often se- 
vere slowing down caused by the substrate to the closest 
layers of the liquid 9-10 . 

In the present work we have considered a much more 
complex system than the bulk. Therefore the fact the 
MCT could be used in this context as an unifying theo- 
retical approach is highly relevant as a guideline for the 
systematic study of the important phenomenology of con- 
fined and interfacial liquids. 
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TABLE I: Parameters of the Lennard- Jones and soft spheres 
potentials as defined in Eq0 In the Table A and B refer 
to the particle of the binary mixture while M refers to the 
confining soft spheres. Values are expressed in Lennard- Jones 
units. 

a e n 
AA 1.0 1.0 1 
BB 0.88 0.5 1 
AB 0.8 1.5 1 
MA 3.0 0.32 
MB 2.94 0.22 
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